Intro

In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.

Load packages & source functions

# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(viridis)
library(sdmTMB)
library(mapplots)
library(sp)
library(raster)
library(ggeffects)
library(modelr)

# Source code for lan_lot to UTM
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/lon_lat_utm.R")

# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/map_plot.R")

# Trim the map plot to better fit stomach data
xmin2 <- 303000
xmax2 <- 916000
xrange <- xmax2 - xmin2

ymin2 <- 5980000
ymax2 <- 6580000
yrange <- ymax2 - ymin2

theme_set(theme_plot())

plot_map_fc2 <- plot_map_fc +
  theme(legend.position = "bottom") +
  xlim(xmin2*1.5, xmax2*0.9) +
  ylim(ymin2*1.025, ymax2*0.98268)

# Continuous colors
options(ggplot2.continuous.colour = "viridis")

Read stomach data

small_cod_stomach <- read_csv("data/clean/small_cod_stomach.csv") %>% 
  drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>% 
  rename(density_cod_small = scod_kg_km2,
         density_cod_large = lcod_kg_km2,
         density_flounder = fle_kg_km2) %>% 
  mutate(fyear = as.factor(year),
         fquarter = as.factor(quarter),
         ices_rect = as.factor(ices_rect),
         oxy_sc = (oxy - mean(oxy)) / sd(oxy),
         temp_sc = (temp - mean(temp)) / sd(temp),
         depth_sc = (depth - mean(depth)) / sd(depth),
         density_saduria_sc = (density_saduria - mean(density_saduria)) / sd(density_saduria),
         density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
         density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
         density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
         haul_id = as.factor(haul_id),
         species2 = "Small cod")
#> Rows: 1153 Columns: 33
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (25): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#> 
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#> 
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#> 
#>         converted 'haul_id' from character to factor (0 new NA)
#> 
#>         new variable 'fyear' (factor) with 6 unique values and 0% NA
#> 
#>         new variable 'fquarter' (factor) with 2 unique values and 0% NA
#> 
#>         new variable 'oxy_sc' (double) with 152 unique values and 0% NA
#> 
#>         new variable 'temp_sc' (double) with 152 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 64 unique values and 0% NA
#> 
#>         new variable 'density_saduria_sc' (double) with 74 unique values and 0% NA
#> 
#>         new variable 'density_cod_small_sc' (double) with 158 unique values and 0% NA
#> 
#>         new variable 'density_cod_large_sc' (double) with 157 unique values and 0% NA
#> 
#>         new variable 'density_flounder_sc' (double) with 159 unique values and 0% NA
#> 
#>         new variable 'species2' (character) with one unique value and 0% NA

large_cod_stomach <- read_csv("data/clean/large_cod_stomach.csv") %>% 
  drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>% 
  rename(density_cod_small = scod_kg_km2,
         density_cod_large = lcod_kg_km2,
         density_flounder = fle_kg_km2) %>% 
  mutate(fyear = as.factor(year),
         fquarter = as.factor(quarter),
         ices_rect = as.factor(ices_rect),
         oxy_sc = (oxy - mean(oxy)) / sd(oxy),
         temp_sc = (temp - mean(temp)) / sd(temp),
         depth_sc = (depth - mean(depth)) / sd(depth),
         density_saduria_sc = (density_saduria - mean(density_saduria)) / sd(density_saduria),
         density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
         density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
         density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
         haul_id = as.factor(haul_id),
         species2 = "Large cod")
#> Rows: 2153 Columns: 33
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (25): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#> 
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#> 
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#> 
#>         converted 'haul_id' from character to factor (0 new NA)
#> 
#>         new variable 'fyear' (factor) with 6 unique values and 0% NA
#> 
#>         new variable 'fquarter' (factor) with 2 unique values and 0% NA
#> 
#>         new variable 'oxy_sc' (double) with 164 unique values and 0% NA
#> 
#>         new variable 'temp_sc' (double) with 164 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 67 unique values and 0% NA
#> 
#>         new variable 'density_saduria_sc' (double) with 76 unique values and 0% NA
#> 
#>         new variable 'density_cod_small_sc' (double) with 158 unique values and 0% NA
#> 
#>         new variable 'density_cod_large_sc' (double) with 171 unique values and 0% NA
#> 
#>         new variable 'density_flounder_sc' (double) with 171 unique values and 0% NA
#> 
#>         new variable 'species2' (character) with one unique value and 0% NA

flounder_stomach <- read_csv("data/clean/flounder_stomach.csv") %>% 
  drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>% 
  rename(density_cod_small = scod_kg_km2,
         density_cod_large = lcod_kg_km2,
         density_flounder = fle_kg_km2) %>% 
  mutate(fyear = as.factor(year),
         fquarter = as.factor(quarter),
         ices_rect = as.factor(ices_rect),
         oxy_sc = (oxy - mean(oxy)) / sd(oxy),
         temp_sc = (temp - mean(temp)) / sd(temp),
         depth_sc = (depth - mean(depth)) / sd(depth),
         density_saduria_sc = (density_saduria - mean(density_saduria)) / sd(density_saduria),
         density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
         density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
         density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
         haul_id = as.factor(haul_id),
         species2 = "Flounder")
#> Rows: 2579 Columns: 33
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (25): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#> 
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#> 
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#> 
#>         converted 'haul_id' from character to factor (0 new NA)
#> 
#>         new variable 'fyear' (factor) with 6 unique values and 0% NA
#> 
#>         new variable 'fquarter' (factor) with 2 unique values and 0% NA
#> 
#>         new variable 'oxy_sc' (double) with 153 unique values and 0% NA
#> 
#>         new variable 'temp_sc' (double) with 153 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 64 unique values and 0% NA
#> 
#>         new variable 'density_saduria_sc' (double) with 77 unique values and 0% NA
#> 
#>         new variable 'density_cod_small_sc' (double) with 148 unique values and 0% NA
#> 
#>         new variable 'density_cod_large_sc' (double) with 160 unique values and 0% NA
#> 
#>         new variable 'density_flounder_sc' (double) with 161 unique values and 0% NA
#> 
#>         new variable 'species2' (character) with one unique value and 0% NA

# Load cache
# qwraps2::lazyload_cache_dir(path = "R/main_analysis/main_diet_analysis_cache/html")

Read the prediction grid (scales variables wrt data!)

# I need predicted cod and flounder densities for this

# pred_grid_1 <- read_csv("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/data/clean/pred_grid_(1_2).csv")
# pred_grid_2 <- read_csv("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/data/clean/pred_grid_(2_2).csv")
# 
# pred_grid <- bind_rows(pred_grid_1, pred_grid_2)

Plot data in space

full_dat <- bind_rows(small_cod_stomach, large_cod_stomach, flounder_stomach) %>% 
  pivot_longer(c(tot_feeding_ratio, saduria_feeding_ratio, benthic_feeding_ratio),
               names_to = "prey_group",
               values_to = "feeding_ratio")
#> pivot_longer: reorganized (tot_feeding_ratio, benthic_feeding_ratio, saduria_feeding_ratio) into (prey_group, feeding_ratio) [was 5885x43, now 17655x42]

## Total feeding ratio
# Q1
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 1 & prey_group == "tot_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt", name = "") +
  labs(title = "Total feeding ratio", subtitle = "Quarter 1") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining


ggsave("figures/supp/data_in_space_q1.pdf", width = 17, height = 17, units = "cm")

# Q4
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 4 & prey_group == "tot_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt", name = "") +
  labs(title = "Total feeding ratio", subtitle = "Quarter 4") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining


ggsave("figures/supp/data_in_space_q4.pdf", width = 17, height = 17, units = "cm")

## Benthic feeding ratio
# Q1
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 1 & prey_group == "benthic_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Benthic feeding ratio", subtitle = "Quarter 1") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining


# Q4
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 4 & prey_group == "benthic_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Benthic feeding ratio", subtitle = "Quarter 4") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining


## Saduria feeding ratio
# Q1
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 1 & prey_group == "saduria_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Saduria feeding ratio", subtitle = "Quarter 1") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining


# Q4
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 4 & prey_group == "saduria_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Saduria feeding ratio", subtitle = "Quarter 4") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining

Plot data (response vs explanatory variables)

ggplot(full_dat, aes(depth, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(oxy, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(temp, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(density_flounder_sc, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(density_cod_large_sc, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(density_cod_small_sc, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


# Saduria
ggplot(full_dat %>% filter(prey_group == "saduria_feeding_ratio"),
       aes(density_saduria_sc, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(~ species2, scales = "free", ncol = 2)
#> filter: removed 11,770 rows (67%), 5,885 rows remaining
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Fit spatiotemporal model to feeding ratio

sdmTMB models with densities as covariates to stomach content data. First make mesh:

# Small cod
small_cod_spde <- make_mesh(small_cod_stomach,
                            c("X", "Y"),
                            cutoff = 15,
                            type = "kmeans",
                            seed = 42)
plot(small_cod_spde)



# Large cod
large_cod_spde <- make_mesh(large_cod_stomach,
                            c("X", "Y"),
                            n_knots = 75,
                            type = "kmeans",
                            seed = 42)
plot(large_cod_spde)



# Flounder
flounder_spde <- make_mesh(flounder_stomach,
                           c("X", "Y"),
                           n_knots = 75,
                           type = "kmeans",
                           seed = 42)
plot(flounder_spde)

# s_cod_tot0 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "off",
#   mesh = small_cod_spde,
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# tidy(s_cod_tot0)
# sanity(s_cod_tot0)
# 
# # Plot residuals in space
# small_cod_stomach$tot_resid0 <- residuals(s_cod_tot0)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid0), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# 
# s_cod_tot <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "off",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# tidy(s_cod_tot)
# sanity(s_cod_tot)
# 
# # Plot residuals in space
# small_cod_stomach$tot_resid <- residuals(s_cod_tot)
# 
# # plot_map_fc2 +
# #   geom_point(data = small_cod_stomach, aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# #              shape = 21, color = "grey50", stroke = 0.2) +
# #   facet_wrap(~ year) +
# #   scale_fill_gradient2()
# #
# # plot_map_fc2 +
# #   geom_point(data = filter(small_cod_stomach, tot_resid < 5), aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# #              shape = 21, color = "grey50", stroke = 0.2) +
# #   facet_wrap(~ year) +
# #   scale_fill_gradient2()
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# # Fit model now with spatial random effect
# s_cod_tot2 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "on",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# print(s_cod_tot2)
# tidy(s_cod_tot2)
# sanity(s_cod_tot2)
# 
# small_cod_stomach$tot_resid2 <- residuals(s_cod_tot2)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid2), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# # Fit model now with spatiotemporal random effect
# s_cod_tot3 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "IID",
#   spatial = "off",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# print(s_cod_tot3)
# tidy(s_cod_tot3)
# sanity(s_cod_tot3)
# 
# small_cod_stomach$tot_resid3 <- residuals(s_cod_tot3)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid3), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# # Fit model now with spatiotemporal AND spatial random effect
# s_cod_tot4 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "IID",
#   spatial = "on",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# print(s_cod_tot4)
# tidy(s_cod_tot4)
# sanity(s_cod_tot4)
# 
# small_cod_stomach$tot_resid4 <- residuals(s_cod_tot4)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid4), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# 
# # Just random effects...
# small_cod_stomach$ices_rect <- as.factor(small_cod_stomach$ices_rect)
# 
# s_cod_tot5 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|ices_rect) + (1|haul_id),
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "off",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# print(s_cod_tot5)
# tidy(s_cod_tot5)
# sanity(s_cod_tot5)
# 
# small_cod_stomach$tot_resid5 <- residuals(s_cod_tot5)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid5), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# 
# # Plot together
# long <- small_cod_stomach %>%
#   pivot_longer(c(tot_resid0, tot_resid, tot_resid2, tot_resid3, tot_resid4, tot_resid5), names_to = "Model", values_to = "Residuals")
# 
# plot_map_fc2 +
#   geom_jitter(data = long %>% filter(Residuals < 4), aes(X*1000, Y*1000, fill = Residuals), size = 3,
#              shape = 21, color = "grey50", stroke = 0.1) +
#   scale_fill_gradient2() +
#   facet_wrap(~Model, ncol = 3)
# 
# plot_map_fc2 +
#   geom_jitter(data = long %>% filter(Residuals < 4), aes(X*1000, Y*1000, fill = Residuals), size = 3,
#              shape = 21, color = "grey50", stroke = 0.1) +
#   scale_fill_gradient2() +
#   facet_grid(year ~ Model)
# 
# 
# # Plot spatial and spatiotemporal random effect
# pred_grid_small <- pred_grid %>%
#   filter(quarter == 4) %>%
#   mutate(X = X/1000,
#          Y = Y/1000) %>%
#   filter(year %in% unique(small_cod_stomach$year)) %>%
#   filter(X < max(small_cod_stomach$X)) %>%
#   filter(X > min(small_cod_stomach$X)) %>%
#   filter(Y < max(small_cod_stomach$Y)) %>%
#   filter(Y > min(small_cod_stomach$Y))
# 
# pred_grid_small <- pred_grid_small %>%
#   mutate(fyear = factor(year),
#          fquarter = factor(quarter),
#          temp_sc = 0,
#          depth_sc = 0,
#          density_cod_small_sc = 0,
#          density_cod_large_sc = 0,
#          density_flounder_sc = 0,
#          oxy_sc = 0)
# 
# spat_pred <- predict(s_cod_tot2, newdata = pred_grid_small)
# 
# plot_map_fc2 +
#   geom_raster(data = spat_pred, aes(X*1000, Y*1000, fill = omega_s)) +
#   scale_fill_gradient2()
# 
# spat_pred2 <- predict(s_cod_tot3, newdata = pred_grid_small)
# 
# plot_map_fc2 +
#   geom_raster(data = spat_pred2, aes(X*1000, Y*1000, fill = epsilon_st)) +
#   scale_fill_gradient2() +
#   facet_wrap(~year)
# 
# AIC(s_cod_tot0, s_cod_tot, s_cod_tot2, s_cod_tot3, s_cod_tot4, s_cod_tot5)

Small cod

# Total
s_cod_tot <- sdmTMB(
  data = small_cod_stomach,
  formula = tot_feeding_ratio ~ 
    0 +
    fyear +
    fquarter +
    temp_sc +
    depth_sc +
    oxy_sc + 
    density_cod_small_sc + 
    density_cod_large_sc + 
    density_flounder_sc +
    (1|haul_id),
  family = delta_lognormal(),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = small_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 6)),
               c(rep(NA, 7), rep(1, 6)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(s_cod_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + 
#>  Formula:     (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#> 
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit') 
#>                      coef.est coef.se
#> fyear2015                1.91    0.64
#> fyear2016                1.53    0.39
#> fyear2017                2.44    0.35
#> fyear2018                1.64    0.30
#> fyear2019                1.47    0.60
#> fyear2020                0.80    0.23
#> fquarter4                0.06    0.52
#> temp_sc                 -0.03    0.16
#> depth_sc                -0.36    0.18
#> oxy_sc                   0.26    0.19
#> density_cod_small_sc     0.01    0.30
#> density_cod_large_sc     0.06    0.29
#> density_flounder_sc      0.20    0.13
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.77
#> 
#> 
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log') 
#>                      coef.est coef.se
#> fyear2015               -4.56    0.34
#> fyear2016               -4.92    0.23
#> fyear2017               -4.34    0.18
#> fyear2018               -5.15    0.18
#> fyear2019               -4.53    0.32
#> fyear2020               -5.00    0.16
#> fquarter4               -0.14    0.27
#> temp_sc                 -0.16    0.10
#> depth_sc                -0.06    0.11
#> oxy_sc                   0.00    0.12
#> density_cod_small_sc    -0.16    0.19
#> density_cod_large_sc     0.11    0.17
#> density_flounder_sc     -0.03    0.08
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.53
#> 
#> Dispersion parameter: 1.26
#> 
#> ML criterion at convergence: -3016.195
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_tot)
#>                    term    estimate std.error
#> 1             fyear2015  1.90743421 0.6411582
#> 2             fyear2016  1.53353316 0.3900081
#> 3             fyear2017  2.43500762 0.3512151
#> 4             fyear2018  1.64370953 0.2964206
#> 5             fyear2019  1.46943783 0.6031040
#> 6             fyear2020  0.80352813 0.2348412
#> 7             fquarter4  0.05884692 0.5173597
#> 8               temp_sc -0.03332226 0.1620087
#> 9              depth_sc -0.36318086 0.1837213
#> 10               oxy_sc  0.26380179 0.1921616
#> 11 density_cod_small_sc  0.01111117 0.3029440
#> 12 density_cod_large_sc  0.06244649 0.2880924
#> 13  density_flounder_sc  0.19657169 0.1297704
sanity(s_cod_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✖ `b_j2` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Benthic
s_cod_ben <- sdmTMB(
  data = small_cod_stomach,
  formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
  family = delta_lognormal(),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = small_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 6)),
               c(rep(NA, 7), rep(1, 6)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(s_cod_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + 
#>  Formula:     (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#> 
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit') 
#>                      coef.est coef.se
#> fyear2015                1.57    0.64
#> fyear2016                1.36    0.39
#> fyear2017                2.15    0.33
#> fyear2018                1.53    0.29
#> fyear2019                0.91    0.60
#> fyear2020                0.68    0.24
#> fquarter4                0.42    0.51
#> temp_sc                 -0.06    0.16
#> depth_sc                -0.38    0.18
#> oxy_sc                   0.39    0.19
#> density_cod_small_sc     0.15    0.30
#> density_cod_large_sc    -0.12    0.27
#> density_flounder_sc      0.22    0.13
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.78
#> 
#> 
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log') 
#>                      coef.est coef.se
#> fyear2015               -4.83    0.34
#> fyear2016               -5.14    0.23
#> fyear2017               -4.54    0.18
#> fyear2018               -5.23    0.18
#> fyear2019               -4.95    0.32
#> fyear2020               -5.21    0.16
#> fquarter4                0.05    0.27
#> temp_sc                 -0.15    0.10
#> depth_sc                -0.05    0.11
#> oxy_sc                   0.15    0.12
#> density_cod_small_sc    -0.03    0.18
#> density_cod_large_sc    -0.04    0.17
#> density_flounder_sc      0.00    0.08
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.52
#> 
#> Dispersion parameter: 1.22
#> 
#> ML criterion at convergence: -2988.002
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_ben)
#>                    term    estimate std.error
#> 1             fyear2015  1.57305072 0.6372890
#> 2             fyear2016  1.35679232 0.3871899
#> 3             fyear2017  2.14593615 0.3327014
#> 4             fyear2018  1.52583561 0.2940587
#> 5             fyear2019  0.91017699 0.5955846
#> 6             fyear2020  0.68246213 0.2358374
#> 7             fquarter4  0.42191915 0.5084646
#> 8               temp_sc -0.05686552 0.1614583
#> 9              depth_sc -0.38373605 0.1808483
#> 10               oxy_sc  0.38603293 0.1892451
#> 11 density_cod_small_sc  0.15307617 0.2985473
#> 12 density_cod_large_sc -0.12246010 0.2670390
#> 13  density_flounder_sc  0.22236500 0.1304145
sanity(s_cod_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j2` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Saduria
s_cod_sad <- sdmTMB(
  data = small_cod_stomach,
  formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*density_saduria_sc + (1|haul_id),
  family = delta_lognormal(),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = small_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 8)),
               c(rep(NA, 7), rep(1, 8)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(s_cod_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc + 
#>  Formula:     depth_sc + density_cod_small_sc + density_cod_large_sc + 
#>  Formula:     density_flounder_sc * density_saduria_sc + (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#> 
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit') 
#>                                        coef.est coef.se
#> fyear2015                                 -1.97    0.83
#> fyear2016                                 -2.33    0.55
#> fyear2017                                 -2.89    0.47
#> fyear2018                                 -2.75    0.48
#> fyear2019                                 -3.49    0.87
#> fyear2020                                 -2.46    0.37
#> fquarter4                                 -0.20    0.68
#> temp_sc                                    0.49    0.23
#> oxy_sc                                     0.57    0.28
#> depth_sc                                   0.41    0.29
#> density_cod_small_sc                       0.11    0.40
#> density_cod_large_sc                       0.02    0.39
#> density_flounder_sc                       -1.12    0.36
#> density_saduria_sc                         0.55    0.23
#> density_flounder_sc:density_saduria_sc     0.59    0.43
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      1.14
#> 
#> 
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log') 
#>                                        coef.est coef.se
#> fyear2015                                 -6.19    0.69
#> fyear2016                                 -5.20    0.47
#> fyear2017                                 -6.21    0.42
#> fyear2018                                 -5.17    0.42
#> fyear2019                                 -5.86    0.78
#> fyear2020                                 -4.87    0.33
#> fquarter4                                  0.17    0.57
#> temp_sc                                   -0.20    0.21
#> oxy_sc                                    -0.16    0.25
#> depth_sc                                   0.07    0.32
#> density_cod_small_sc                      -0.18    0.50
#> density_cod_large_sc                       0.27    0.66
#> density_flounder_sc                        0.20    0.36
#> density_saduria_sc                         0.27    0.21
#> density_flounder_sc:density_saduria_sc     0.30    0.42
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.14
#> 
#> Dispersion parameter: 1.30
#> 
#> ML criterion at convergence: -271.004
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_sad)
#>                                      term    estimate std.error
#> 1                               fyear2015 -1.96975855 0.8299213
#> 2                               fyear2016 -2.33073173 0.5472003
#> 3                               fyear2017 -2.88815853 0.4729243
#> 4                               fyear2018 -2.75489048 0.4813353
#> 5                               fyear2019 -3.49083692 0.8668283
#> 6                               fyear2020 -2.45756405 0.3749937
#> 7                               fquarter4 -0.19681017 0.6823595
#> 8                                 temp_sc  0.49499626 0.2292552
#> 9                                  oxy_sc  0.56663816 0.2773261
#> 10                               depth_sc  0.40580797 0.2873256
#> 11                   density_cod_small_sc  0.10921056 0.4002075
#> 12                   density_cod_large_sc  0.01993056 0.3870146
#> 13                    density_flounder_sc -1.11983774 0.3647768
#> 14                     density_saduria_sc  0.55029773 0.2316145
#> 15 density_flounder_sc:density_saduria_sc  0.58664443 0.4258442
sanity(s_cod_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Test residuals for single model, tweedie vs lognormal
# s_cod_tot2 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + 
#     density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
#   family = tweedie(),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "off",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# # Delta-Lognormal
# small_cod_stomach$res <- residuals(s_cod_tot)
# 
# p1 <- plot_map_fc2 +
#   geom_point(data = small_cod_stomach, aes(X*1000, Y*1000, fill = res), size = 2, shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2() + 
#   ggtitle("Delta-Lognormal")
# 
# mcmc_res <- residuals(s_cod_tot, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
# p1b <- ggplot(as.data.frame(mcmc_res), aes(sample = mcmc_res)) + stat_qq() + stat_qq_line() + theme(aspect.ratio = 1)
# qqnorm(mcmc_res);qqline(mcmc_res)
# 
# # Tweedie
# small_cod_stomach$res2 <- residuals(s_cod_tot2)
# 
# p2 <- plot_map_fc2 +
#   geom_point(data = small_cod_stomach, aes(X*1000, Y*1000, fill = res2), size = 2, shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2() + 
#   ggtitle("Tweedie")
# 
# mcmc_res2 <- residuals(s_cod_tot2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200) %>% as.data.frame()
# p2b <- ggplot(mcmc_res2, aes(sample = V1)) + stat_qq() + stat_qq_line() + theme(aspect.ratio = 1)
# 
# # Plot!
# (p1 + p2 + p1b + p2b) + plot_layout(heights = c(2, 1))

Large cod

# Total
l_cod_tot <- sdmTMB(
  data = large_cod_stomach,
  formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
  family = delta_lognormal(),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = large_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 6)),
               c(rep(NA, 7), rep(1, 6)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(l_cod_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + 
#>  Formula:     (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#> 
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit') 
#>                      coef.est coef.se
#> fyear2015                1.38    0.35
#> fyear2016                1.49    0.23
#> fyear2017                1.96    0.20
#> fyear2018                1.34    0.16
#> fyear2019                1.06    0.35
#> fyear2020                1.22    0.17
#> fquarter4                0.38    0.31
#> temp_sc                 -0.46    0.11
#> depth_sc                -0.11    0.13
#> oxy_sc                   0.05    0.13
#> density_cod_small_sc     0.76    0.30
#> density_cod_large_sc    -0.60    0.30
#> density_flounder_sc      0.00    0.08
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.44
#> 
#> 
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log') 
#>                      coef.est coef.se
#> fyear2015               -3.08    0.33
#> fyear2016               -3.80    0.21
#> fyear2017               -3.44    0.17
#> fyear2018               -3.58    0.16
#> fyear2019               -3.71    0.32
#> fyear2020               -3.74    0.17
#> fquarter4               -0.79    0.27
#> temp_sc                  0.04    0.10
#> depth_sc                 0.03    0.13
#> oxy_sc                  -0.05    0.12
#> density_cod_small_sc     0.20    0.28
#> density_cod_large_sc    -0.05    0.28
#> density_flounder_sc     -0.09    0.08
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.55
#> 
#> Dispersion parameter: 1.76
#> 
#> ML criterion at convergence: -4845.201
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_tot)
#>                    term     estimate  std.error
#> 1             fyear2015  1.376267080 0.35248671
#> 2             fyear2016  1.489466984 0.23264291
#> 3             fyear2017  1.955719828 0.20279435
#> 4             fyear2018  1.342176159 0.16264333
#> 5             fyear2019  1.063579950 0.35163094
#> 6             fyear2020  1.224150901 0.17030026
#> 7             fquarter4  0.382795458 0.31010761
#> 8               temp_sc -0.459569859 0.11054982
#> 9              depth_sc -0.109393256 0.13361705
#> 10               oxy_sc  0.054455940 0.12583284
#> 11 density_cod_small_sc  0.757360174 0.30113449
#> 12 density_cod_large_sc -0.597163459 0.29919650
#> 13  density_flounder_sc  0.002566405 0.07747939
sanity(l_cod_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Benthic
l_cod_ben <- sdmTMB(
  data = large_cod_stomach,
  formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
  family = delta_lognormal(),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = large_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 6)),
               c(rep(NA, 7), rep(1, 6)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(l_cod_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + 
#>  Formula:     (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#> 
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit') 
#>                      coef.est coef.se
#> fyear2015               -0.18    0.36
#> fyear2016                0.15    0.22
#> fyear2017                0.41    0.18
#> fyear2018               -0.09    0.16
#> fyear2019               -0.78    0.34
#> fyear2020                0.02    0.17
#> fquarter4                1.26    0.30
#> temp_sc                 -0.38    0.11
#> depth_sc                -0.10    0.14
#> oxy_sc                   0.28    0.13
#> density_cod_small_sc     0.25    0.30
#> density_cod_large_sc    -0.21    0.29
#> density_flounder_sc      0.05    0.08
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.58
#> 
#> 
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log') 
#>                      coef.est coef.se
#> fyear2015               -5.50    0.31
#> fyear2016               -5.71    0.20
#> fyear2017               -5.56    0.17
#> fyear2018               -5.63    0.16
#> fyear2019               -6.35    0.30
#> fyear2020               -5.72    0.17
#> fquarter4                0.39    0.26
#> temp_sc                 -0.12    0.10
#> depth_sc                -0.04    0.13
#> oxy_sc                   0.08    0.12
#> density_cod_small_sc    -0.41    0.27
#> density_cod_large_sc     0.52    0.27
#> density_flounder_sc     -0.09    0.07
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.49
#> 
#> Dispersion parameter: 1.55
#> 
#> ML criterion at convergence: -4771.129
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_ben)
#>                    term    estimate  std.error
#> 1             fyear2015 -0.17736775 0.35712898
#> 2             fyear2016  0.15284441 0.21906882
#> 3             fyear2017  0.41442854 0.18416163
#> 4             fyear2018 -0.09346653 0.15820153
#> 5             fyear2019 -0.78207442 0.34474478
#> 6             fyear2020  0.01629230 0.17278533
#> 7             fquarter4  1.26366931 0.29920972
#> 8               temp_sc -0.38340307 0.10941454
#> 9              depth_sc -0.09719543 0.13672989
#> 10               oxy_sc  0.28408572 0.12649178
#> 11 density_cod_small_sc  0.25424771 0.29604417
#> 12 density_cod_large_sc -0.20528910 0.29446028
#> 13  density_flounder_sc  0.04657151 0.08037281
sanity(l_cod_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Saduria
l_cod_sad <- sdmTMB(
  data = large_cod_stomach,
  formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*density_saduria_sc + (1|haul_id),
  family = delta_lognormal(),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = large_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 8)),
               c(rep(NA, 7), rep(1, 8)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(l_cod_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc + 
#>  Formula:     depth_sc + density_cod_small_sc + density_cod_large_sc + 
#>  Formula:     density_flounder_sc * density_saduria_sc + (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#> 
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit') 
#>                                        coef.est coef.se
#> fyear2015                                 -1.80    0.83
#> fyear2016                                 -3.39    0.54
#> fyear2017                                 -3.64    0.51
#> fyear2018                                 -3.65    0.44
#> fyear2019                                 -4.15    0.92
#> fyear2020                                 -2.82    0.43
#> fquarter4                                 -0.05    0.71
#> temp_sc                                    0.00    0.25
#> oxy_sc                                     0.39    0.31
#> depth_sc                                   0.23    0.36
#> density_cod_small_sc                       0.33    0.56
#> density_cod_large_sc                      -0.38    0.60
#> density_flounder_sc                       -0.83    0.31
#> density_saduria_sc                         0.40    0.23
#> density_flounder_sc:density_saduria_sc     0.40    0.38
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      1.53
#> 
#> 
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log') 
#>                                        coef.est coef.se
#> fyear2015                                 -5.58    0.54
#> fyear2016                                 -6.02    0.35
#> fyear2017                                 -5.74    0.36
#> fyear2018                                 -5.95    0.29
#> fyear2019                                 -6.02    0.71
#> fyear2020                                 -6.07    0.32
#> fquarter4                                 -0.78    0.53
#> temp_sc                                    0.53    0.18
#> oxy_sc                                     0.33    0.24
#> depth_sc                                   0.35    0.31
#> density_cod_small_sc                      -0.16    0.50
#> density_cod_large_sc                       0.28    0.77
#> density_flounder_sc                       -0.12    0.22
#> density_saduria_sc                         0.16    0.18
#> density_flounder_sc:density_saduria_sc    -0.12    0.34
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.39
#> 
#> Dispersion parameter: 1.29
#> 
#> ML criterion at convergence: -498.453
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_sad)
#>                                      term     estimate std.error
#> 1                               fyear2015 -1.799469107 0.8330093
#> 2                               fyear2016 -3.385412445 0.5368501
#> 3                               fyear2017 -3.640778227 0.5114008
#> 4                               fyear2018 -3.650167817 0.4403155
#> 5                               fyear2019 -4.147620658 0.9194658
#> 6                               fyear2020 -2.821250957 0.4257319
#> 7                               fquarter4 -0.052218529 0.7129595
#> 8                                 temp_sc  0.002607801 0.2518307
#> 9                                  oxy_sc  0.389346476 0.3110741
#> 10                               depth_sc  0.227847260 0.3615738
#> 11                   density_cod_small_sc  0.328026403 0.5551598
#> 12                   density_cod_large_sc -0.382032051 0.5989242
#> 13                    density_flounder_sc -0.834920812 0.3136644
#> 14                     density_saduria_sc  0.397871281 0.2278797
#> 15 density_flounder_sc:density_saduria_sc  0.396319690 0.3848843
sanity(l_cod_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

Flounder

# Total
fle_tot <- sdmTMB(
  data = flounder_stomach,
  formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
  family = delta_lognormal(),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = flounder_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 6)),
               c(rep(NA, 7), rep(1, 6)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(fle_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + 
#>  Formula:     (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#> 
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit') 
#>                      coef.est coef.se
#> fyear2015                0.53    0.53
#> fyear2016                0.30    0.39
#> fyear2017                0.36    0.25
#> fyear2018                0.12    0.25
#> fyear2019               -0.03    0.47
#> fyear2020                1.42    0.23
#> fquarter4                0.86    0.38
#> temp_sc                 -0.10    0.14
#> depth_sc                -0.20    0.17
#> oxy_sc                   0.33    0.19
#> density_cod_small_sc     0.63    0.37
#> density_cod_large_sc    -0.58    0.37
#> density_flounder_sc      0.03    0.12
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      1.04
#> 
#> 
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log') 
#>                      coef.est coef.se
#> fyear2015               -4.37    0.24
#> fyear2016               -4.40    0.18
#> fyear2017               -3.98    0.13
#> fyear2018               -4.04    0.13
#> fyear2019               -4.06    0.22
#> fyear2020               -3.64    0.10
#> fquarter4                0.37    0.18
#> temp_sc                  0.18    0.07
#> depth_sc                 0.00    0.09
#> oxy_sc                   0.13    0.10
#> density_cod_small_sc    -0.19    0.19
#> density_cod_large_sc     0.21    0.19
#> density_flounder_sc     -0.07    0.05
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.38
#> 
#> Dispersion parameter: 1.26
#> 
#> ML criterion at convergence: -3530.324
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_tot)
#>                    term    estimate std.error
#> 1             fyear2015  0.52611178 0.5277953
#> 2             fyear2016  0.29633519 0.3855055
#> 3             fyear2017  0.36424865 0.2499490
#> 4             fyear2018  0.11872342 0.2529731
#> 5             fyear2019 -0.03306316 0.4684866
#> 6             fyear2020  1.42321774 0.2260591
#> 7             fquarter4  0.86332478 0.3823155
#> 8               temp_sc -0.10027050 0.1439928
#> 9              depth_sc -0.20081704 0.1677329
#> 10               oxy_sc  0.32846294 0.1869042
#> 11 density_cod_small_sc  0.62940403 0.3726341
#> 12 density_cod_large_sc -0.57761745 0.3693939
#> 13  density_flounder_sc  0.03169148 0.1150707
sanity(fle_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✖ `b_j2` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Benthic
fle_ben <- sdmTMB(
  data = flounder_stomach,
  formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
  family = delta_lognormal(),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = flounder_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 6)),
               c(rep(NA, 7), rep(1, 6)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(fle_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + 
#>  Formula:     (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#> 
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit') 
#>                      coef.est coef.se
#> fyear2015                0.51    0.54
#> fyear2016                0.27    0.39
#> fyear2017                0.32    0.25
#> fyear2018                0.14    0.26
#> fyear2019               -0.28    0.48
#> fyear2020                1.33    0.23
#> fquarter4                0.86    0.39
#> temp_sc                 -0.05    0.15
#> depth_sc                -0.22    0.17
#> oxy_sc                   0.35    0.19
#> density_cod_small_sc     0.66    0.38
#> density_cod_large_sc    -0.62    0.37
#> density_flounder_sc      0.05    0.12
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      1.06
#> 
#> 
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log') 
#>                      coef.est coef.se
#> fyear2015               -4.39    0.23
#> fyear2016               -4.43    0.17
#> fyear2017               -3.98    0.13
#> fyear2018               -4.01    0.13
#> fyear2019               -4.14    0.22
#> fyear2020               -3.65    0.10
#> fquarter4                0.37    0.18
#> temp_sc                  0.19    0.07
#> depth_sc                -0.07    0.09
#> oxy_sc                   0.09    0.10
#> density_cod_small_sc    -0.21    0.18
#> density_cod_large_sc     0.23    0.18
#> density_flounder_sc     -0.05    0.05
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.35
#> 
#> Dispersion parameter: 1.28
#> 
#> ML criterion at convergence: -3476.879
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_ben)
#>                    term    estimate std.error
#> 1             fyear2015  0.50568356 0.5372639
#> 2             fyear2016  0.27473734 0.3927066
#> 3             fyear2017  0.31836479 0.2545055
#> 4             fyear2018  0.13679684 0.2572819
#> 5             fyear2019 -0.27996769 0.4768295
#> 6             fyear2020  1.32914073 0.2284407
#> 7             fquarter4  0.86119225 0.3891565
#> 8               temp_sc -0.04676748 0.1462127
#> 9              depth_sc -0.21504409 0.1708393
#> 10               oxy_sc  0.35217170 0.1903804
#> 11 density_cod_small_sc  0.65658216 0.3776414
#> 12 density_cod_large_sc -0.61939795 0.3746410
#> 13  density_flounder_sc  0.05410617 0.1171719
sanity(fle_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Saduria
fle_sad <- sdmTMB(
  data = flounder_stomach,
  formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc + depth_sc + density_flounder_sc + density_cod_large_sc + density_cod_small_sc*density_saduria_sc + (1|haul_id),
  family = delta_lognormal(),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = flounder_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 8)),
               c(rep(NA, 7), rep(1, 8)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(fle_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc + 
#>  Formula:     depth_sc + density_flounder_sc + density_cod_large_sc + density_cod_small_sc * 
#>  Formula:     density_saduria_sc + (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#> 
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit') 
#>                                         coef.est coef.se
#> fyear2015                                  -1.08    0.85
#> fyear2016                                  -0.50    0.65
#> fyear2017                                  -2.43    0.45
#> fyear2018                                  -3.65    0.49
#> fyear2019                                  -2.69    0.82
#> fyear2020                                  -2.07    0.40
#> fquarter4                                  -0.31    0.64
#> temp_sc                                    -0.20    0.23
#> oxy_sc                                      0.26    0.30
#> depth_sc                                    0.07    0.27
#> density_flounder_sc                        -0.42    0.26
#> density_cod_large_sc                        0.57    0.54
#> density_cod_small_sc                       -0.82    0.58
#> density_saduria_sc                          0.68    0.26
#> density_cod_small_sc:density_saduria_sc    -0.11    0.52
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      1.75
#> 
#> 
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log') 
#>                                         coef.est coef.se
#> fyear2015                                  -3.92    0.41
#> fyear2016                                  -4.00    0.30
#> fyear2017                                  -5.18    0.24
#> fyear2018                                  -4.79    0.29
#> fyear2019                                  -4.09    0.43
#> fyear2020                                  -4.31    0.22
#> fquarter4                                  -1.06    0.31
#> temp_sc                                    -0.15    0.13
#> oxy_sc                                     -0.57    0.19
#> depth_sc                                   -0.49    0.20
#> density_flounder_sc                        -0.39    0.15
#> density_cod_large_sc                        0.24    0.36
#> density_cod_small_sc                       -0.21    0.43
#> density_saduria_sc                          0.29    0.15
#> density_cod_small_sc:density_saduria_sc     0.16    0.33
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.53
#> 
#> Dispersion parameter: 1.27
#> 
#> ML criterion at convergence: -1139.288
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_sad)
#>                                       term    estimate std.error
#> 1                                fyear2015 -1.08052722 0.8501067
#> 2                                fyear2016 -0.50181784 0.6485335
#> 3                                fyear2017 -2.43386977 0.4479369
#> 4                                fyear2018 -3.64549688 0.4871173
#> 5                                fyear2019 -2.68587056 0.8176042
#> 6                                fyear2020 -2.07202390 0.3970650
#> 7                                fquarter4 -0.31306727 0.6417368
#> 8                                  temp_sc -0.19566564 0.2337761
#> 9                                   oxy_sc  0.26092780 0.3040494
#> 10                                depth_sc  0.06517983 0.2721824
#> 11                     density_flounder_sc -0.41502450 0.2621885
#> 12                    density_cod_large_sc  0.56748090 0.5371357
#> 13                    density_cod_small_sc -0.81664754 0.5787022
#> 14                      density_saduria_sc  0.68338619 0.2593120
#> 15 density_cod_small_sc:density_saduria_sc -0.10661946 0.5234617
sanity(fle_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
mcmc_res_fle_sad <- residuals(fle_sad,
                              type = "mle-mcmc",
                              mcmc_iter = 201,
                              mcmc_warmup = 200)
#> Constructing atomic log_dbinom_robust
#> Constructing atomic invpd
#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001647 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 16.47 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 14.2302 seconds (Warm-up)
#> Chain 1:                0.026116 seconds (Sampling)
#> Chain 1:                14.2563 seconds (Total)
#> Chain 1:

qqnorm(mcmc_res_fle_sad); qqline(mcmc_res_fle_sad)

Plot coefficients

#### Total prey
# small cod total prey
small_cod_tot_pars <- bind_rows(tidy(s_cod_tot, effects = "fixed", conf.int = TRUE, model = 1) %>% 
                                  mutate(model = "Binomial"),
                                tidy(s_cod_tot, effects = "fixed", conf.int = TRUE, model = 2) %>% 
                                  mutate(model = "Lognormal")) %>%
  mutate(species = "Small cod",
         response = "Total prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# large cod total prey
large_cod_tot_pars <- bind_rows(tidy(l_cod_tot, effects = "fixed", conf.int = TRUE, model = 1) %>% 
                                  mutate(model = "Binomial"),
                                tidy(l_cod_tot, effects = "fixed", conf.int = TRUE, model = 2) %>% 
                                  mutate(model = "Lognormal")) %>% 
  mutate(species = "Large cod",
         response = "Total prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# flounder total prey
flounder_tot_pars <- bind_rows(tidy(fle_tot, effects = "fixed", conf.int = TRUE, model = 1) %>% 
                                  mutate(model = "Binomial"),
                                tidy(fle_tot, effects = "fixed", conf.int = TRUE, model = 2) %>% 
                                  mutate(model = "Lognormal")) %>% 
  mutate(species = "Flounder",
         response = "Total prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

#### Benthic prey
# small cod benthic prey
small_cod_ben_pars <- bind_rows(tidy(s_cod_ben, effects = "fixed", conf.int = TRUE, model = 1) %>% 
                                  mutate(model = "Binomial"),
                                tidy(s_cod_ben, effects = "fixed", conf.int = TRUE, model = 2) %>% 
                                  mutate(model = "Lognormal")) %>% 
  mutate(species = "Small cod",
         response = "Benthic prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# large cod benthic prey
large_cod_ben_pars <- bind_rows(tidy(l_cod_ben, effects = "fixed", conf.int = TRUE, model = 1) %>% 
                                  mutate(model = "Binomial"),
                                tidy(l_cod_ben, effects = "fixed", conf.int = TRUE, model = 2) %>% 
                                  mutate(model = "Lognormal")) %>% 
  mutate(species = "Large cod",
         response = "Benthic prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# flounder benthic prey
flounder_ben_pars <- bind_rows(tidy(fle_ben, effects = "fixed", conf.int = TRUE, model = 1) %>% 
                                  mutate(model = "Binomial"),
                                tidy(fle_ben, effects = "fixed", conf.int = TRUE, model = 2) %>% 
                                  mutate(model = "Lognormal")) %>% 
  mutate(species = "Flounder",
         response = "Benthic prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA


#### Saduria
# small cod saduria
small_cod_sad_pars <- bind_rows(tidy(s_cod_sad, effects = "fixed", conf.int = TRUE, model = 1) %>% 
                                  mutate(model = "Binomial"),
                                tidy(s_cod_sad, effects = "fixed", conf.int = TRUE, model = 2) %>% 
                                  mutate(model = "Lognormal")) %>% 
  mutate(species = "Small cod",
         response = "Saduria")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# large cod saduria
large_cod_sad_pars <- bind_rows(tidy(l_cod_sad, effects = "fixed", conf.int = TRUE, model = 1) %>% 
                                  mutate(model = "Binomial"),
                                tidy(l_cod_sad, effects = "fixed", conf.int = TRUE, model = 2) %>% 
                                  mutate(model = "Lognormal")) %>% 
  mutate(species = "Large cod",
         response = "Saduria")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# flounder saduria
flounder_sad_pars <- bind_rows(tidy(fle_sad, effects = "fixed", conf.int = TRUE, model = 1) %>% 
                                  mutate(model = "Binomial"),
                                tidy(fle_sad, effects = "fixed", conf.int = TRUE, model = 2) %>% 
                                  mutate(model = "Lognormal")) %>% 
  mutate(species = "Flounder",
         response = "Saduria")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

coefs <- bind_rows(small_cod_tot_pars,
                   large_cod_tot_pars,
                   flounder_tot_pars,
                      
                   small_cod_ben_pars,
                   large_cod_ben_pars,
                   flounder_ben_pars,
                   
                   small_cod_sad_pars,
                   large_cod_sad_pars,
                   flounder_sad_pars)

coefs <- coefs %>% 
  mutate(term2 = term,
         term2 = ifelse(term == "temp_sc", "Temperature", term2),
         term2 = ifelse(term == "oxy_sc", "Oxygen", term2),
         term2 = ifelse(term == "depth_sc", "Depth", term2),
         term2 = ifelse(term == "density_saduria_sc", "Saduria", term2),
         term2 = ifelse(term == "density_cod_large_sc", "Cod (> 29)", term2),
         term2 = ifelse(term == "density_cod_small_sc", "Cod (< 29)", term2),
         term2 = ifelse(term == "density_cod_small_sc:density_saduria_sc", "Cod (> 29) x Saduria", term2),
         term2 = ifelse(term == "density_flounder_sc:density_saduria_sc", "Flounder x Saduria", term2),
         term2 = ifelse(term == "density_cod_small_sc:density_saduria_sc", "Small cod x Saduria", term2),
         term2 = ifelse(term == "density_flounder_sc", "Flounder", term2))
#> mutate: new variable 'term2' (character) with 16 unique values and 0% NA

coefs <- coefs %>% 
  mutate(sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
         sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))
#> mutate: new variable 'sig' (character) with 2 unique values and 0% NA

# Fixed continuous effects:
# coefs %>% 
#   filter(!grepl('year', term)) %>% 
#   filter(!grepl('quarter', term)) %>% 
#   ggplot(aes(term2, estimate, color = species)) +
#   geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
#   geom_point(position = position_dodge(width = 0.4)) +
#   geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
#                 position = position_dodge(width = 0.4), alpha = 0.5) +
#   scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
#   facet_grid(model~response) +
#   labs(x = "Predictor", y = "Standardized coefficient") +
#   theme(legend.position = "bottom",
#         axis.text.x = element_text(angle = 90),
#         aspect.ratio = 1) +
#   NULL 

#ggsave("figures/fr_coefs.pdf", width = 17, height = 17, units = "cm")

coefs %>% 
  filter(!grepl('year', term)) %>% 
  filter(!grepl('quarter', term)) %>% 
  ggplot(aes(term2, estimate, color = model, fill = model, shape = sig)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
                position = position_dodge(width = 0.4), alpha = 0.5) +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  scale_shape_manual(values = c(1, 21)) +
  guides(fill = "none") +
  facet_grid(response ~ species) +
  labs(x = "", y = "Standardized coefficient", color = "Model", shape = "95% CI crossing 0") +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90, vjust = 0.4)) +
  NULL 
#> filter: removed 108 rows (44%), 138 rows remaining
#> filter: removed 18 rows (13%), 120 rows remaining


ggsave("figures/fr_coefs.pdf", width = 17, height = 17, units = "cm")

# Year effects
coefs %>% 
  filter(term %in% c("fyear2015", "fyear2016", "fyear2017", "fyear2018", "fyear2019", "fyear2020", "quarter")) %>% 
  mutate(term2 = ifelse(term == "fyear2015", "2015", term2),
         term2 = ifelse(term == "fyear2016", "2016", term2),
         term2 = ifelse(term == "fyear2017", "2017", term2),
         term2 = ifelse(term == "fyear2018", "2018", term2),
         term2 = ifelse(term == "fyear2019", "2019", term2),
         term2 = ifelse(term == "fyear2020", "2020", term2),
         term2 = ifelse(term == "quarter", "Quarter", term2)) %>% 
  filter(!term2 == "Quarter") %>% 
  ggplot(aes(term2, estimate, color = species)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  facet_grid(model~response) +
  labs(x = "Predictor", y = "Standardized coefficient") +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90),
        aspect.ratio = 1.5) +
  NULL 
#> filter: removed 138 rows (56%), 108 rows remaining
#> mutate: changed 108 values (100%) of 'term2' (0 new NA)
#> filter: no rows removed


ggsave("figures/supp/fr_coefs_yr.pdf", width = 17, height = 17, units = "cm")
# fle_sad <- sdmTMB(
#   data = flounder_stomach,
#   formula = saduria_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "off",
#   mesh = flounder_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )

# visreg_delta(s_cod_sad, xvar = "density_flounder_sc", model = 1, gg = TRUE, by = "oxy_sc")
# visreg_delta(s_cod_sad, xvar = "density_flounder_sc", model = 2, gg = TRUE)
# visreg_delta(s_cod_sad, xvar = "density_flounder_sc", model = 1, gg = TRUE, by = "oxy_sc")


hist(fle_sad$data$density_cod_small_sc)


fle1 <- visreg_delta(fle_sad,
                     xvar = "density_cod_small_sc",
                     model = 1,
                     plot = FALSE,
                     by = "density_saduria_sc",
                     breaks = c(-0.5, 0.5)
                     )

fle2 <- visreg_delta(fle_sad,
                     xvar = "density_cod_small_sc",
                     model = 2,
                     plot = FALSE,
                     by = "density_saduria_sc",
                     breaks = c(-0.5, 0.5))

scod1 <- visreg_delta(s_cod_sad,
                      xvar = "density_flounder_sc",
                      model = 1,
                      plot = FALSE,
                      by = "density_saduria_sc",
                      breaks = c(-0.5, 0.5))

scod2 <- visreg_delta(s_cod_sad,
                      xvar = "density_flounder_sc",
                      model = 2,
                      plot = FALSE,
                      by = "density_saduria_sc",
                      breaks = c(-0.5, 0.5))

lcod1 <- visreg_delta(l_cod_sad,
                      xvar = "density_flounder_sc",
                      model = 1,
                      plot = FALSE,
                      by = "density_saduria_sc",
                      breaks = c(-0.5, 0.5))

lcod2 <- visreg_delta(l_cod_sad,
                      xvar = "density_flounder_sc",
                      model = 2,
                      plot = FALSE,
                      by = "density_saduria_sc",
                      breaks = c(-0.5, 0.5))

visreg_fit <- bind_rows(fle1$fit %>% mutate(model = "fle1"),
                        fle2$fit %>% mutate(model = "fle2"),
                        scod1$fit %>% mutate(model = "scod1"),
                        scod2$fit %>% mutate(model = "scod2"),
                        lcod1$fit %>% mutate(model = "lcod1"),
                        lcod2$fit %>% mutate(model = "lcod2")
                        )
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA

visreg_res <- bind_rows(fle1$res %>% mutate(model = "fle1"),
                        fle2$res %>% mutate(model = "fle2"),
                        scod1$res %>% mutate(model = "scod1"),
                        scod2$res %>% mutate(model = "scod2"),
                        lcod1$res %>% mutate(model = "lcod1"),
                        lcod2$res %>% mutate(model = "lcod2")
                        )
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA

visreg_fit <- visreg_fit %>% 
  mutate(species = ifelse(model %in% c("fle1", "fle2"), "Flounder", NA),
         species = ifelse(model %in% c("scod1", "scod2"), "Small cod", species),
         species = ifelse(model %in% c("lcod1", "lcod2"), "Large cod", species)) %>% 
  mutate(model_comp = ifelse(model %in% c("fle1", "scod1", "lcod1"), "Binomial", "LogNormal"))
#> mutate: new variable 'species' (character) with 3 unique values and 0% NA
#> mutate: new variable 'model_comp' (character) with 2 unique values and 0% NA

visreg_res <- visreg_res %>% 
  mutate(species = ifelse(model %in% c("fle1", "fle2"), "Flounder", NA),
         species = ifelse(model %in% c("scod1", "scod2"), "Small cod", species),
         species = ifelse(model %in% c("lcod1", "lcod2"), "Large cod", species)) %>% 
  mutate(model_comp = ifelse(model %in% c("fle1", "scod1", "lcod1"), "Binomial", "LogNormal"))
#> mutate: new variable 'species' (character) with 3 unique values and 0% NA
#> mutate: new variable 'model_comp' (character) with 2 unique values and 0% NA

pal <- brewer.pal(n = 9, name = "Set1")[c(2, 9)]

# Flounder
p1 <- ggplot(visreg_fit %>% filter(species == "Flounder"),
             aes(x = density_cod_small_sc, y = visregFit, fill = factor(density_saduria_sc), color = factor(density_saduria_sc))) +
  facet_grid(model_comp ~ species) +
  geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2, color = NA) +
  geom_line() +
  labs(x = "Small cod density (scaled)",
       color = "Saduria density (scaled)") +
  guides(fill = "none") +
  geom_point(aes(y = visregRes), data = visreg_res %>% filter(species == "Flounder"), size = 1, alpha = 0.5) + 
  # scale_color_brewer(palette = "Set1") +
  # scale_fill_brewer(palette = "Set1") +
  # scale_color_brewer(palette = "Accent") +
  # scale_fill_brewer(palette = "Accent") +
  scale_color_manual(values = pal) +
  scale_fill_manual(values = pal) +
  theme(aspect.ratio = 1.1) +
  NULL
#> filter: removed 808 rows (67%), 404 rows remaining
#> filter: removed 3,658 rows (54%), 3,174 rows remaining

# Small cod
p2 <- ggplot(visreg_fit %>% filter(species == "Small cod"),
             aes(x = density_flounder_sc, y = visregFit, fill = factor(density_saduria_sc), color = factor(density_saduria_sc))) +
  #ggh4x::facet_grid2(model_comp ~ species, scales = "free_y", independent = "y") +
  facet_grid(model_comp ~ species) +
  geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2, color = NA) +
  geom_line() +
  labs(x = "Flounder density (scaled)",
       color = "Saduria density (scaled)") +
  guides(fill = "none") +
  geom_point(aes(y = visregRes), data = visreg_res %>% filter(species == "Small cod"), size = 1, alpha = 0.2) + 
  scale_color_manual(values = pal) +
  scale_fill_manual(values = pal) +
  theme(aspect.ratio = 1.1) +
  NULL
#> filter: removed 808 rows (67%), 404 rows remaining
#> filter: removed 5,536 rows (81%), 1,296 rows remaining

# Large cod
p3 <- ggplot(visreg_fit %>% filter(species == "Large cod"),
             aes(x = density_flounder_sc, y = visregFit, fill = factor(density_saduria_sc), color = factor(density_saduria_sc))) +
  #ggh4x::facet_grid2(model_comp ~ species, scales = "free_y", independent = "y") +
  #ggh4x::facet_grid2(model_comp ~ species, scales = "free_y", independent = "y") +
  facet_grid(model_comp ~ species) +
  geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2, color = NA) +
  geom_line() +
  labs(x = "Flounder density (scaled)",
       color = "Saduria density (scaled)") +
  guides(fill = "none") +
  geom_point(aes(y = visregRes), data = visreg_res %>% filter(species == "Large cod"), size = 1, alpha = 0.2) + 
  scale_color_manual(values = pal) +
  scale_fill_manual(values = pal) +
  theme(aspect.ratio = 1.1) +
  NULL
#> filter: removed 808 rows (67%), 404 rows remaining
#> filter: removed 4,470 rows (65%), 2,362 rows remaining

(p1 | p2 | p3) + plot_layout(guides = "collect")

  
# visreg_fit %>%
#   #filter(species == "Large cod") %>%
#   group_by(species, model) %>%
#   summarise(max_covar = max(density_flounder_sc),
#             min_covar = min(density_flounder_sc))

# Option 2
# ggplot(visreg_fit, aes(x = density_flounder_sc, y = visregFit, fill = factor(density_saduria_sc), color = factor(density_saduria_sc))) +
#   facet_grid2(model_comp ~ species) +
#   geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2, color = NA) +
#   geom_line() +
#   labs(x = "Flounder density (scaled)",
#        color = "Saduria density") +
#   guides(fill = "none") +
#   geom_point(aes(y = visregRes), data = visreg_res, size = 1, alpha = 0.2) + 
#   scale_color_brewer(palette ="Dark2") +
#   scale_fill_brewer(palette = "Dark2") +
#   theme(aspect.ratio = 1.1) +
#   NULL

ggsave("figures/conditional_saduria.pdf", width = 17, height = 17, units = "cm")